Full characterization of the fractional Poisson process 
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The fractional Poisson process (FPP) is a counting process with independent and identically 
distributed inter-event times following the Mittag-Lefner distribution. This process is very useful in 
several fields of applied and theoretical physics including models for anomalous diffusion. Contrary 
to the well-known Poisson process, the fractional Poisson process does not have stationary and 
independent increments. It is not a Levy process and it is not a Markov process. In this letter, we 
present formulae for its finite-dimensional distribution functions, fully characterizing the process. 
These exact analytical results are compared to Monte Carlo simulations. 



From a loose mathematical point of view, counting pro- 
cesses N(A) are stochastic processes that count the ran- 
dom number of points in a set A. They are used in many 
fields of physics and other applied sciences. In this let- 
ter, we will consider one-dimensional real sets with the 
physical meaning of time intervals. The points will be in- 
coming events whose duration is much smaller than the 
inter-event or inter-arrival waiting time. For instance, 
counts from a Geiger-Muller counter can be described in 
this way. The number of counts, N(At), in a given time 
interval At is known to follow the Poisson distribution 

P(N(At) = n) = exp(-AAt) (AA ,^" , (1) 

TV. 

where A is the constant rate of arrival of ionizing par- 
ticles. Together with the assumption of independent 
and stationary increments, Eq. ([!]) is sufficient to de- 
fine the homogeneous Poisson process. Curiously, one 
of the first occurrences of this process in the scientific 
literature was connected to the number of casualties by 
horse kicks in the Prussian army cavalry pQ. The Pois- 
son process is strictly related to the exponential dis- 
tribution. The inter-arrival times Ti identically follow 
the exponential distribution and are independent ran- 
dom variables. This means that the Poisson process is 
a prototypical renewal process. A justification for the 
ubiquity of the Poisson process has to do with its rela- 
tionship with the binomial distribution. Suppose that 
the time interval of interest (t, t + At) is divided into 
n equally spaced sub-intervals. Further assume that a 
counting event appears in such a sub-interval with prob- 
ability p and does not appear with probability 1 — p. 
Then, ¥(N(At) = k) = Bin(fc;p, n) is a binomial distri- 
bution of parameters p and n and the expected number 
of events in the time interval is given by E[iV(Ai)] = np. 
If this expected number is kept constant for n — > oo, 
the binomial distribution converges to the Poisson dis- 
tribution of parameter A = E[N(At)}/ At, while, in the 



meantime, p — > 0. However, it can be shown that many 
counting processes with non-stationary increments con- 
verge to the Poisson process after a transient period. It 
is sufficient to require that they are renewal process (i.e. 
they have independent and identically distributed (iid) 
inter-arrival times) and that E(Tj) < oo. In other words, 
many counting processes with non-independent and non- 
stationary increments behave as the Poisson process if 
observed long after the transient period. 

In recent times, it has been shown that heavy-tailed 
distributed inter-arrival times (for which E(rj) = oo) do 
play a role in many phenomena such as blinking nano- 
dots [2 |3] , human dynamics [U [5] and the related inter- 
trade times in financial markets [51 [7] . 

Among the counting processes with non-stationary in- 
crements, the so-called fractional Poisson process [H], 
Np(t), is particularly important because it is the thinning 
limit of counting processes related to renewal processes 
with power-law distributed inter-arrival times [9l 110] . 
Moreover, it can be used to approximate anomalous dif- 
fusion ruled by space-time fractional diffusion equations 
[H [TTHTtj] . It is a straightforward generalization of the 
Poisson process defined as follows. Let {n}^ 1 be a se- 
quence of independent and identically distributed posi- 
tive random variables with the meaning of inter-arrival 
times and let their common cumulative distribution func- 
tion (cdf) be 

F T (t)=F(r<t) = l-E (-t% (2) 

where Ep{—t^) is the one-parameter Mittag-Leffler func- 
tion, Ep(z), defined in the complex plane as 

oo n 

evaluated in the point z = —t^ and with the prescription 
< ft < 1. In equation T(-) is Euler's Gamma 



2 



function. The sequence of the epochs, {T n }^ 
by the sums of the inter-arrival times 



=1) is given 



(4) 



The epochs represent the times in which events arrive 
or occur. Let f T (t) = dF T (t)/dt denote the probability 
density function (pdf) of the inter-arrival times, then the 
probability density function of the n-th epoch is simply 
given by the n-fold convolution of f T (t), written as f* n (t). 
In Ref. [TO], it is shown that 



f Tn (t) = f* T n (t) = P 



t nfi~ 



(n - 1) 



] E 



(5) 



where (— i* 3 ) is the n-th derivative of Ep{z) evalu- 
ated in z = —t 13 . The counting process Np(t) counts the 
number of epochs (events) up to time t, assuming that 
T = is an epoch as well, or, in other words, that the 
process begins from a renewal point. This assumption 
will be used all over this paper. Np(t) is given by 



Np(t) = max{n : T„ < t}. 



(6) 



In Ref. [S], the fractional Poisson distribution is derived 
and it is given by 



_P=0.5, t 1= 5 
P=0-8, t.=5 
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P(Np(t) = n) 



{-If). 



(7) 



Eq. Q coincides with the Poisson distribution of param- 
eter A = 1 for j3 = 1. In principle, equations ^ and 
([7]) can be directly used to derive the fractional Poisson 
distribution, but convergence of the series is slow. Fortu- 
nately, in a recent paper, Beghin and Orsingher proved 
that 



Figure 1: P(Np(Ti) = m) as function of ni for three different 
values of (3. The crosses are estimations obtained from 10 s 
Monte Carlo samples and the lines are given to guide the eye. 



Fs„(t;u) 



exp(—u)u r ' 
(n-l)l 



exp(— u)u" 



du. 



(8) 



where Fg^ (t; u) is the cdf of a stable random variable 
S/3(u, 7, S) with index /3, skewness parameter v = 1, scale 
parameter 7 = (ucostt/3/2) 1 /^ and location (5 = [T7] , 
The integral in equation ^ can be evaluated numerically 
and Fig. [T]l shows P(Np(t) = n) for three different val- 
ues of j3. The Monte Carlo simulation of the fractional 
Poisson process is based on the algorithm presented in 
equation (20) of Ref. [H]. 

As a consequence of Kolmogorov's extension theorem, 
in order to fully characterize the stochastic process Np (t) , 
one has to derive its finite dimensional distributions. A 
further requirement on the process' paths uniquely deter- 
mines the process, namely that they are right-continuous 
step functions with left limits |18j . The finite-dimensional 



Y„, 



Figure 2: (Color online) Pictorial illustration of the random 
variables used in the text. The light blue dots represent the 
observation points t\, t-z and t^. The red squares are the 
epochs To = 0, Ti, . . . ,Ts- The conditional residual life-time 
is the time elapsed between ti and the next epoch T ni +i. It 
depends on previous values of ni, this is the number of events 
between and ti, with the event at t = To = not considered. 
Here, we have m = 1, 712 = 2 and n-i — 4. All the equations 
in this paper can be derived by analyzing this figure. 
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distributions are the multivariate probability distribution 
functions V(Np(ti) = n x , Np(t 2 ) = n 2 , . . . , Np(t k ) = n k ) 
with t\ < t 2 < ... < t k and n\ < n 2 < •■■ < n k . 
We have already given the formula for the one-point 
functions in Eq. Q. The general finite dimensional 
distribution can be computed observing that the event 
{Np(ti) = ni,Np(t 2 ) = n 2 , . . .,Np(t k ) = n k } is equiva- 
lent to {0 < T ni < ti,T„ 1+ i > ti,ti < T n2 < t 2 ,T„ 2+ i > 
t 2 , . . .,ifc_i < T nk < t k ,T nk+1 > t k }. Therefore, we find 

F{Np(t 1 )=n 1 ,Np(t 2 ) = n 2 ,...,Np(t k )=n k ) = 
P(0 < T ni < ti,T ni+ i > t\,t\ < T„ 2 < t 2 ,T n2+ i > t 2 , 
■ ■ ■ ,tk-x < T nh < t k ,T nk+ i > t k ) = 

d Ul / T * ni («i) / du 2 f T (u 2 ) 



t-i—ui 



t\—Ui—U2 



dtt3 /;(n 2 -«i-l)( n3 ) / du 4 f T (u 4 ) 

J t 2 — U 1 —U2—U 3 



du 2k ^f< n "- n "-^{u 2k ^ 

2k-l 



th-l— 5Di=l 2 u 



l-F T [t k -Y J 



Hi 



i=l 



(9) 



For instance, the two point function is given by 



F(Np(ti)=n 1 ,Np(t 2 ) = n 2 ) = 

< T ni < ti,T, n+ i > ti,ti < T n2 < t 2 ,T, l2+ i > t 2 ) 



dmf^im) / du 2 f T (u 2 ) 

Jtl—Ul 
*2 — Ml — U 2 

d U3 /;(«-»i-i)( U 3) 
[1 - F T (ta - ui - ua - us)] . (10) 

Let us focus on the two-point case for the sake of 
illustration. As Np(t) is a counting process, one has 
¥{N p {tx) = n x ,Np{t 2 ) = n 2 ) = P(Np(h) = n u Np{t 2 ) - 
Np(t±) = n 2 —ni) and, as a consequence of the definition 
of conditional probability 

V(Np(h) = m,Np{t 2 ) - Np(t 1 )=n 2 -n 1 ) = 
F{N p (t 2 ) - Np(h) =n 2 - rulNpih) = m)x 

xP(JV j8 (ti)=ni). (11) 

For (3 — 1, when the fractional Poisson process coincides 
with the standard Poisson process, the increments are iid 
random variables and one has 

P(JVi(* a ) - Ntih) =n 2 - ni|JV x (ti) - n x ) = 
P(JVi(t 2 )-JVi(ti) =n 2 -m) = 

expH^-tO) ^ 2 lj . (12) 

(n 2 -ni)! 



On the contrary, for < f3 < 1, the increment Np(t 2 ) — 
Np(t\) and Np(t\) are not independent. Note that 
Np(t\) can be seen as an increment as Np(0) — by 
definition. However from Eq. (11 1, the conditional prob- 
ability of having n 2 — ri\ epochs in the interval {t\,t 2 ) 
conditional on the observation of n\ epochs in the interval 
(0, t\) can be written as a ratio of two finite dimensional 
distribution: 



F(Np(t 2 ) - Np(h) = na - n^N^h) = m) = 
F(Np(h) = n u Np(t 2 ) = n 2 ) 



F(Np(t 1 )=n 1 ) 



(13) 



This probability can be evaluated by means of an alter- 
native method, more appealing for a direct and practical 
understanding of the dependence structure. Let 



del' 



Y ni = [T ni+1 -t 1 \Np(t 1 )=n 1 



(14) 



denote the residual lifetime at time t\ (that is the time to 
the next epoch or renewal) conditional on Np(t±) = ?ii . 
With reference to Fig. [2j one can see that the conditional 
probability F(Np(t 2 ) - Np(ti) =n 2 - ni\Np(ti) = ni) is 
given by the following convolution integral for n 2 — ri\ > 1 

P(Np(t 2 ) - Np{h) =n 2 - n^Npih) = ni) = 



F(Np(t 2 -h-y)= n 2 - m - l)f Yni (y) dy, 



(15) 



where fy n (t) is the pdf of Y ni . In the case n 2 — = 0, 
one has 

P(Np(t 2 ) - N p (ti) = 0|JV^(ti) = m) = 1 - F Yni (i a - h) 

(16) 

where Fy ni (y) is the cdf of Y ni . The distribution of the 
conditional residual lifetime Y ni can be evaluated in sev- 
eral ways. For instance, one can notice that it can be 
decomposed as follows 



Y n 1 T n 1 _|_ i -\- U n 1 



where U ni is defined as 



U ni = [T ni \Np{t 1 )=n 1 ], 



(17) 



(18) 



and is the position of the last epoch before t\ conditional 
on Np(t\) = ni, and 



T ni + 1 — [ T n 1 + 1 — tl\T ni+ i > fi] 



(19) 



is the difference between t„ 1+ i and t\ conditional on 
T ni +i > t\. The pdf of U ni is given by the following 



chain of equalities 



where we defined 



f Vni (t)dt =P(T ni g dt\N p (ti) = m) 

=P(T ni € dt|T ni < ti,T ni +r ni+1 > h] 
=P(T m G dt|T ni < ii,T Wl+ i > *i - T ni ; 
, P(T ttl edt)/~_ t P(r Wl+1 efa) 

P(r ni < ti,T„ 1+ i > ti — T ni ) 

* ^(t)[l-f r (tl-t)]dt 



5^n,...,n* d = P^+i - tk\Np(ti) = ni, . . . , A^(tfc) = rife]. 

(27) 

Again, we can use a decomposition of y n i,...,n fc 



ft dufr'Hil-Frih-u)} 7 



(20) 



where we used the independence between T ni and t„ 1+ i 
(*) and / Tbj (Z) = f*^ (x) (*). The pdf of f ni+1 is 



f? ni+ x(t\U ni )dt =P(t„ 1+ i - ti G dt|T ni+ i > t x ) 

= P(r ni +i €j* + fr) 
P(r„ 1+ i > *i-l7 ni ) 



(21) 



From Eq. (17), one can write that 



ff ni+1 (t-u\u)f Uni {u)du (22) 



and this equation leads to 

rti 



ft dufi^MMt + h-u) 

ft du/r x («)[i--FV(*i-u)] 



that, together with Eq. Q, gives us the probability of 
the conditional increments in Eq. (15). Notice that, for 
rii = 0, one has f*°(u) = S(u) and Eq. (23) reduces to 



the familiar equation for the residual life-time pdf in the 
absence of previous renewals 



JVo(t) 



fr(t + h) 



(24) 



l-F T (h) 

This method can be applied to the general multidimen- 



sional case. As in Eq. (11) we can write 



P(iV/j(*i) = n u ...,Ni3(tk) = n k ,N p (t k+1 ) = n h+1 ) = 
®>(Np(t k+l ) - Npfa) = n k+i ~ n k \ 
Npfa) = ni, . . . , Np(t k ) = n k )x 

xF(Np(t 1 )=n 1> ...,Np(t k )=n k ) (25) 

and the predictive probabilities can be evaluated as 

¥(N p (t k+1 ) - Np{t k ) = n fc+ i - n k \ . . . 

\Np{t 1 ) = n 1 ,...,Np{t k )=n k ) = 



W{N p {t k+1 -tk-y) = n k+ i - n k - l)x 



x fY ni ...., nk (y)dy, (26) 



(28) 



where 



U nk = [TnjNpfa) =n u .. .,N p {t k ) = n k ], (29) 



and 



tlef 



Tn k + 1 - [r nh +i - t k T„.+l > t k 



(30) 



The difference with the two-point case is that U ni — 
[TnMfa) = m] = [YZ± 1 n\N fl (t 1 ) = ni] must be re- 
placed by 



^ n\Np(t k ) = n k 



i=n fc -i+l 



(31) 

The time between t k -\ and the next renewal epoch 
is Y ni ,...,n k „ l and it is independent from E^n+i^- 
Therefore, the convolution 

q( ni ,...,n k ;t) = f Yni „^ * £(»*-»»-i-i)( t ) (32) 



(23) replaces f* ni (t) in Eq. ((201). This leads to 



g(ni,...,n fc ;f + f fc _i)[l-F r fa - f)] 
/ t ^ g(ni, . . . ,n k ;u + t fc _i)[l - F T (t k - u)]di 



(33) 



On the other hand, ff n +1 (t) has the same functional 



form as /f n +1 (t) given in Eq. (21) with L/ nfc replacing 
J7 ni . Therefore, Y nit __ nh has the following pdf 

J t *_ i dixg(ni,...,nfc;u + tfc_i)/ r (f + tA - ti) 



■^fe-i dM( ?(«i, ■ ■ • + *fe-i)[i — ^V(* + *fc — «)] 



(34) 



carries the 



In practice, the random variable Y ni t1 
memory of the observations made at times t±, . . . , t k -\] 
the knowledge of fy ni „ j allows the computation 
°f „ ' an dj via Eqs. (25) and (26), the k + 1- 



dimensional distribution can be derived as well. 

Figs. P and g] compare the theoretical results of 
Eqs. (|20E p3l and pi} with those of a Monte Carlo 



simulation based on the algorithm presented in equation 
(20) of Ref. PJ]. 
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Figure 3: (Color online) Pdf of the random variable U ni as given in Eq. (201 (solid black lines) compared to Monte Carlo 
simulations (colored step lines) for three values of /3 and two different values of t\. 10 7 different paths were simulated for each 
value of /3 and the bin width is 0.05. Time is in arbitrary units. 



G 






Figure 4: (Color online) Pdf of the random variable Y ni as given in Eqs. (231 and (241 (solid black lines) compared to Monte 
Carlo simulations (colored step lines) for three values of (3 and two different values of ti . 10 7 different paths were simulated for 
each value of /3 and the bin width is 0.01. Time is in arbitrary units. 
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